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Abstract 

The article addresses the important for homeland security possibility of robust 
detection of geometrically small, low emission sources on a significantly stronger back- 
ground. A technique of detecting such sources using Compton type cameras is devel- 
oped and studied analytically and numerically, which is shown to have high sensitivity 
and specificity and also allows to assign confidence probabilities of the detection. The 
2D case is considered in detail. 

1 Introduction 

One of the missions of the Department of Homeland Security is to prevent smuggling 
of weapon-grade nuclear materials (e.g., |25|). It is expected that such materials, 
unlike those needed for a "dirty bomb", will have low emission rates and will be well 
shielded, so that very few gamma photons or neutrons would escape, and even less 
would be detected non-scattered (ballistic). An additional hurdle for the detection of 
illicit nuclear substances is the strong natural radiation background. If the particles 
radiated by the source were physically distinct (e.g., in terms of their energies) from the 
ones prevalent in the background, then the detection would be easier. We thus assume 
that the particles from the source are identical to those coming from the background 
in terms of their nature, energies, etc., so the discrimination using such parameters 
is impossible. This leads to the situation when the signal to noise ratio (SNR), i.e. 
the count of the detected non-scattered particles emitted by the source versus the 
total number of detected particles, can be as low as 10~^, if not lower. It is rather 
clear that if only the locations of the hits are detected, with no information about 
the directions of the incoming particles available, there is no chance to detect the 
presence of such a weak source. Thus, one needs to employ detectors that can provide 
some directional information. For instance, the 7-cameras most commonly used in 
SPECT (Single Photon Emission Computed Tomography) medical imaging |4,27 , are 
mechanically collimated, and thus they "count" the particles coming from a narrow 
cone of directions and discard the rest. However, mechanical collimation dramatically 
reduces the particle count, and thus, in the case of an extremely low emission source, 
can essentially eliminate the useful signal (this difficulty, although in much less extreme 
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form, arises also in SPECT). Typically, only one out of thousands emitted photons is 
detected by a collimated camera. When dealing with low emission sources, one wants 
to capture as many particles as possible, and thus mechanical collimation is unsuitable 
for the problem at hand. 

Another option is to use the so called Compton 7-cameras (see Section [5]), which 
do not discard any incoming particles by collimation. The price one pays for this is 
that the directional information provided by such a camera is more limited. Namely, 
one obtains just a hollow (i.e., surface rather than solid) cone of possible directions 
rather than a single direction. Still, Compton cameras are good candidates for the 
applications we have in mind. While Compton cameras are available for detecting 
7-photons, neutron detectors that provide a similar "incoming cone" information are 
being currently developed as well |25ll40|. The methods we will discuss do not depend 
upon a specific kind of particles; we will thus call all detectors that can provide the 
cone information Compton type detectors. 

Even if some directional information is available, the task of detecting extremely 
low emission coming from a geometrically large source in the presence of a significant 
radiation background would be very hard, if not impossible. So, another important 
condition that we impose, besides availability of directional information, is the small 
size of the possible source, which is usually a safe assumption in the applications we 
are interested in. 

Simulations of radiation from a small high-enriched uranium source placed in a cargo 
container suggest that only about 0.1% of the signal received by detectors might be due 
to the ballistic (non-scattered) particles emitted from the source. The remaining 99.9% 
of detected particles come either from extraneous (natural) sources, or from scattered 
source particles from the source |5|. 

It is natural to try to use the rather standard SPECT techniques in the present 
situation, albeit the chances of success (and even the applicability of the tomographic 
models) are questionable. This circle of issues is discussed analytically and tested 
numerically in Section [2j where the conclusion is made that backprojection technique 
might be the best bet here, while the usual filtration parts in SPECT reconstructions do 
not do any good. This conclusion gets its foundation in the probabilistic discussions of 
Section |3] and then in numerical examples of Section [4| However, to check the viability 
of the approach, all these considerations dealt with the simplest, collimated detectors. 
Section [5] is devoted to the case of Compton type cameras. Here tomographic and 
backprojection techniques are considered. The analysis of the 2D Compton camera 
case is provided in Section [6j The corresponding numerical tests are conducted and 
discussed in Section [7| The overall conclusion is that the suggested backprojection 
Compton type cameras techniques allow a robust detection of presence of geometrically 
small low emission sources with extremely low SNR, with confidence levels attached. 

Like in the case of SPECT |4|, one can try to use statistical methods. E.g., the 
recent paper |46| describes a Bayesian approach to the same problem. 

The paper ends with the sections devoted to remarks and conclusions, acknowledg- 
ments, and bibliography. 
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2 Source detection using tomographic techniques 



We will briefly present here some relevant mathematical formulas and conclusions from 
the SPECT version of emission tomography, which is the closest to the problem of our 
interest. One can find more details in g[[2l|[22l[27[[2| . 

Let f{x) be the intensity distribution of the sources of particles of certain type 
(7-photons, neutrons, etc.) inside an object (e.g., cargo container, or a truck). Let 
also fi{x) be the attenuation coefficient (usually not known the applications we have 
in mind). Then, in the approximation of sufficiently high emission rate, the radia- 
tion transport equation implies that the particle count per unit of time at a detector 
collimated in the direction L is 

T^f{L)^ jj{x)e-^L.^(y)dydx. (1) 

Here Lx is the segment of the line L between the emission point x and the detector 
and dy denotes the standard linear measure on L. The operator is said to be the 
attenuated Radon (or X-ray) transform (with attenuation //(x)) of the function 




Figure 1: The set-up of the SPECT emission tomography 



We will fix an origin O and assign the normal coordinates (c<j, s) to any line L. Here 
uj — (cos 0, sin 6) is a unit vector orthogonal to L and s is the (signed) distance from the 
origin to L (see Fig. [T]). The equation of this line is x • uo — s and the line itself will be 
sometimes denoted as L^^^^^y The value of ([T]) provides the mathematical expectation 
of the counts at a detector per unit of time, which is not a good approximation of the 
measured data, if the emission rate is low. Let us still accept for the time being this 
integral-geometric formulation and see where it leads us. 

The attenuated backprojection operator with attenuation v{x^uj) > (which 
may or may not be related to fi{x)) acting on data ^(cj, s) is 

[Tfg){x)^l g{uj,x • u)v{x,uj)duj (2) 

^ ^ J\uj\=l 

In particular, ii v — 1 (which we will assume in our experiments), this is the standard 
backprojection operator [20||27||28] . 
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Let H be the Hilbert transform acting in variable s on functions g{uj^ s)\ 

TT J S- 



Here v.p. indicates that the integral is considered in the principal value meaning. 
The following result holds: 

Theorem 1. (see details in [221) Let f{x) be a function (our distribution of radiating 
sources) and fi be sufficiently smootl^ then the operator 



T*hIt, (3) 



applied to f preserves all singularities of f and does not create new ones. In particular, 
if f has a jump across an interface, the location of this jump is preserved in the image 

Notice that Tf H-^T^f is not a reconstruction of / from its measured projections 
T^/ (which would be impossible anyway with the low signal strength and the attenu- 
ation ji being unknown), but it gives a correct reconstruction of all singularities (e.g., 
interfaces) of /. 

Since we assume that our source is very small geometrically, it is rather singular, 
and thus one might expect that any guess for the unknown attenuation (e.g., that there 
is none) should reconstruct the location of the source correctly. 

The operator 

T*hI (4) 

is exactly what is applied to the data in our experiments described below. We thus 

will not worry about the presence of an unknown attenuation. 

The operator (W) is often called a filtered backprojection (FBP), where H-^ 
41. ' — ' 

and are the filtration and backprojection parts, correspondingly. 

Remark 2. • Since we are looking for a rather singular object, it is reasonable to 
recall that there is a way in tomography to emphasize singularities of the function 
to be reconstructed. This is what the so called local tomography does (e.g., 
1^ ~2^). For instance, the singularities will look brighter if one replaces the 
filtration part H in with a stronger filter 

This option will also be explored. 

• In the detection problem we are discussing, different types of particles (e.g., 7- 
photons of different energies) will be present, and each of them could have different 
distribution fj of sources, as well as different attenuation distribution fij in the 
cargo. Thus, what we get if we lump all these particle counts together (which 
seems to be a strange thing to do), is the sum 



^The detection procedures we will eventually come to, will not require smoothness of attenuation. 
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Can this complicate things? We claim that it should not. Indeed^ all these func- 
tions have the same boundary of the source as their interface. 

One can derive now from the above theorem the following 

Corollary 3. Applying the same reconstruction operator T^H-^ to the whole 
sum, one recovers the correct source interface. 

It is still not clear whether such a procedure will work in our problem, since there 
are at least two potentially serious obstacles: 

• The (attenuated) X-ray transform model does not apply when significant statis- 
tical noise is present, in particular for a very weak source. 

• This model does not take into account presence of a very large (from standard 
tomographic point of view, enormous) background noise. 

Still, we will try to apply this approach and see where it can lead us. We show 
below some sample results of numerical experimentation with the 2D and 3D X-ray 
transforms, as well as 3D Radon transform. 

2.1 2D X-ray reconstructions 

Although the problem is three-dimensional, for the first trial we consider the 2D prob- 
lem, which is more straightforward to handle numerically. 

In the first reconstruction we have modeled random 7-photon emission from a 
shielded ball of radius 4cm of HEU235 for realistic values of the emission and at- 
tenuation rates. The ball was placed inside in the centeij^of a much (about 100 times) 
larger "cargo". A uniform and isotropic random background was also added. The 
model included 64 detectors, each of which could detect 64 directional sectors (about 
2.8 degrees each). In this reconstruction we assumed the SNR to be 1%, i.e. 99% 
of detected hits were generated by background, or scattered source particles. Then, 
the described above tomographic X-ray transform inversion (see ([s]) with v — 1) was 
applied to the resulting matrix of counts. The pictures below show the results of the 
reconstructions in 44 x 44 pixels (both the gray scale density and surface plots are 
presented). The result, shown in Fig. [2| shows clear detection of the source. 

We then try to lower the SNR towards our benchmark value 0.1%, and the success 
seems to evaporate. E.g., Fig. [3] shows the failed attempt of reconstruction with the 
SNR being 0.4% and with only about a hundred ballistic particles detected from the 
source. 

Now let us consider a fixed signal-to- noise ratio (0.1%) and varying total number of 
detected particles. We will first describe the set up for the numerical results presented 
in rest of this subsection. 

Four direction sensitive detector arrays, each consisting of 100 detectors (bins), were 
placed along the sides of the square [—1,1]^. Random background radiation was created 
by simulating particles propagating along straight lines in such a way that these lines 
were uniformly randomly distributed in the square. More precisely, to each particle 
we assigned a direction of propagation by choosing the normal coordinates (6^, s) of 
a line such that 9 and s were uniformly distributed in (— 7r/2,7r/2) and [— \/2, \/2], 
correspondingly. Lines which pass through the imaged square necessarily intersect two 

^The results for non-central locations of the source, as the further examples will show, are similar. 
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Figure 2: Reconstruction with SNR 1% (density and surface plots). 




Figure 3: Unsuccessful FBP reconstruction with SNR 0.4% and about 100 ballistic source 
particles detected (density and surface plots). 



detector arrays, and one of the two intersection points was randomly chosen as the site 
of detection of the particle. Finally, for each such particle we recorded the bin on the 
detector array and the exact incoming direction a (see Fig [4]). We also simulated an 
isotropic point source which emitted a number of particles roughly equal to 0.1% of 
the background. In order to use the filtered backprojection (FBP) inversion formula 
Q we first computed the normal coordinates of the (straight-line) trajectories of the 
particles: 

= a — 7r/2, s — [xi^yi) • (cos0, sin0), 

where (x^, yi) is the center of the bin on the detector array. After that, the coordinates 
6 and s were put into bins of size 0.02. A two dimensional data array was constructed 
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in which every element equaled the number of particles with trajectories in the cor- 
responding 0, s bin. Different tomographic inversions were applied to such data. The 
criterion for successful detection was whether the highest pick of the reconstructed 
image occurred in a small neighbourhood of the source location. 




Figure 4: A particle propagates along a random line with normal coordinates (6^, s) and 
collides with a detector array at the point (x, y). We record the bin with center (x^, yi) and 
the direction a = 9 + 7t/2. 



In our first experiment, we simulate about 300,000 background particles and ad- 
ditional 300 particles coming from a source located at (0.401, —0.133). In most cases, 
filtered backprojection reconstruction Q would fail to detect the source: only 7 out of 
20 experiments were successful. However, increasing the total number of detected par- 
ticles, while keeping the SNR at 0.1%, improved the results significantly. Setting the 
background to 400, 000 and then to 500, 000 particles resulted in 13 and 18 detections 
out of 20 experiments, correspondingly. 

As we have already mentioned, the local tomography (which uses a stronger 
filter, see ([s])) emphasizes the singularities of an image, and thus one could expect it 
to detect geometrically small sources better. This, however, turns out not to be the 
case. In our experiments local tomography always showed inferior results to filtered 
backprojection, often failing to detect a source when an FBP detection from the same 
data was successful. A stronger high-pass filter just increases the contribution from the 
noise, so the reconstructed images, in fact, worsened. This suggests that maybe one 
should try to eliminate the filtration step altogether and relay upon the backprojection 
alone. 

Indeed, our results using backprojection with no filtration were very encouraging. 
Backprojection of the same data used for FBP and local reconstructions resulted in 
16, 18 and 20 out of 20 successful detections for 300, 000, 400, 000 and 500, 000 detected 
background particles, correspondingly. The results of these experiments are summa- 
rized in Table |2.1[ Fig [5] shows an example of an unsuccessful FBP reconstruction 
and the corresponding successful backprojection reconstruction from data generated 
by 499, 765 background particles and 500 source paricles. 
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method \ bkgd particles 


300, 000 


400, 000 


500, 000 


Backprojection 
FBP 

Local tomography 


16/20 
7/20 
6/20 


18/20 
13/20 
9/20 


20/20 
18/20 
13/20 



Table 1: Number of successful detections of a source located at (0.401,-0.133), with a 
fixed SNR 0.1% and varying number of detected particles, using diflFerent reconstruction 
methods. For each level of background particles 20 sets of random data were generated. In 
each case backprojection, FBP and local tomography were used to detect the source. In all 
experiments, if a local tomography detection was successful, so was the corresponding FBP 
detection, and if there was an FBP detection, backprojection would also find the source. 




Figure 5: Reconstructions from x-ray data generated by 499, 765 random background parti- 
cles and 500 particles emitted by an isotropic source located at (0.401, —0.133). The black 
lines indicate the location of detector arrays. The correct source location is encircled. Top 
row: FBP reconstruction (left), only peaks deviating more than 3.2 standard deviations 
from the mean (right). A number of peaks in the image are higher than the peak at the 
source. Bottom row: backprojection reconstruction (left), only peaks deviating more than 
3.5 standard deviations from the mean (right). The correct location of the source is found. 
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2.2 3D X-ray and Radon reconstructions 



While in 2D there is a single transform that can be called interchangeably Radon 
transform or X-ray transform, these notions part in higher dimensions. In 3D, there is 
the X-ray transform, which integrates a function f{x) over all straight lines 
in R^, as well as the Radon transform, which integrates functions over all 2D 

planes. Both are of interest in tomography, the first arising in X-ray CT, SPECT, 
and PET scanners and the second in MRI imaging. In the situation we consider in this 
text, X-ray transform represents, the same way as it does in 2D, the case of collimated 
detectors, while the Radon transform will have some relation to the Compton camera 
case considered in Section [Sj Indeed, in both 3D Radon and Compton cases, the data 
represents surface integrals (over planes and cones respectively). 



2.2.1 3D X-ray reconstructions 

The same issues that we have discussed in 2D case remain: the X-ray transform model 
is valid only for strong sources (or long observation time) and it does not address 
the background noise. However, in 3D there is another feature that one has to take 
into account. Namely, while in 2D the set of lines is also 2-dimensional (the lines 
can be parametrized by their two normal coordinates), in 3D the space of lines is 4- 
dimensional. Thus, the X-ray data in 3D is redundant. In tomography, this redundancy 
is usually dealt with by selecting a 3-dimensional sub-set of data (the rest might not 
even been collected). For instance, one can use only rays parallel to a fixed plane 
and do the reconstruction slice-by-slice, or one can pick only lines intersecting a given 
curve (often a helix). In the detection situation we are facing, this is not an option, 
since any attempt to reduce dimension of the data will essentially remove all signal 
and leave only noise in its wake. Thus, we need to use all, overdetermined data. This 
requires corresponding analytic formulas (which are easy to obtain and thus will not 
be introduced) and much heavier calculations. 

Without further ado, we present an example (Fig. [g]) of a successful FBP recon- 
struction, where 10^ ballistic particles from the source and 1.35 x 10^ background ones 
were detected. 

Further experimenting with examples leads to the same conclusion as in 2D: back- 
projection is superior to FBP in this situation, while local reconstruction leads to 
deterioration. 



2.2.2 3D Radon reconstructions 

In this section, we present some Radon transform reconstructions in 3D. The reader 
should recall that the 3D Radon transform integrates functions over affine 2D planes, 
rather than lines. This is a preliminary test of the principle before doing 3D recon- 
structions for Compton cameras, where surface integrals (over cones) are also involved. 

In the figures [t]- [9] below, there were 10^ background particles detected, while the 
number of ballistic particles from the source was varying, thus changing the SNR and 
the level of detectability. Density and surface plots of 2D sections through the source 
locations are shown. One sees that with the SNR of 0.5% the source can be detected, 
while at 0.1% and lower the detection deteriorates. 



Let us see now whether backprojection alone does a better job. Figures [T0| and 11 



show the backprojection reconstruction from the same data as in Fig. [8] and |9| One 
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Figure 6: Here one sees a "transparent" 3D visualization, as well as pieces of ID and 2D 
sections passing through the source. 




Figure 7: FBP reconstruction from 3D Radon data. 5000 ballistic source particles. SNR 
0.5%. 



clearly sees that the backprojection alone detects the source when the FBP method 
fails to do so. 

Again, the examples above confirm the general feasibility of the approach and 
superiority of the backprojection. 

2.3 Discussion of the tomographic reconstructions 

The experiments described in the previous sections beg several questions: 
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Figure 8: FBP reconstruction from 3D Radon data. 1000 ballistic source particles. SNR 
0.1%. 




Figure 9: FBP reconstruction from 3D Radon data. 700 ballistic source particles. SNR 
0.07%. 



1. Why do tomographic methods still work to some degree, although the assump- 
tions that lead to the X-ray transform model are not satisfied (i.e., the source is 
too weak) and a strong random background is present? Indeed, the use of stan- 
dard tomography assumes that the data represents line integrals of the source 
distribution function. This is the same as to say that we measure the expectation 
of the number of hits, per unit time, from particles moving along a line. This 
assumption is reasonable when a substantial number of source particles is de- 
tected, i.e. when the source is strong, or the observation time is long. However, it 
fails in the case of low emission sources that cannot be observed for a really long 
time. Hence, the X-ray/Radon transform type models and techniques appear to 
be inappropriate. 

2. Why does strengthening the filter (i.e., local tomography), which should increase 
detectability of singularities, makes the reconstruction worse? 

3. Why do tomographic techniques fail to reliably detect sources at the levels of 
SNR and total number of particles for which the probabilistic consideration of 
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Figure 10: BP reconstruction from 3D Radon data. 1000 ballistic source particles, SNR 
0.1%. Top: Cross section through source. Bottom: Isosurfaces at 50% (blue), 60% (green) 
and 70% (red) of peak value. 



the next section suggests that sources should be detectable? 

As we might have guessed already, the answer to the first question is that in fact, 
tomographic inversions do not work. Or, to put it differently, only a part of the 
algorithm (backprojection) works, while its other part (filtration) only hurts. 

The answer to the second question is not hard to guess. Indeed, high-pass filters 
required in X-ray/Radon transform inversions amplify noise. This is a very significant 
factor, as the noise constitutes 99.9% of our data. Using local tomography can only 
aggravate this difficulty, and we have seen that it indeed does. 

The answer to the third question is probably the same as the previous one: high- 
pass filters might be the culprits. 

These considerations suggest to try to eliminate high-pass filtering completely, and 
thus to use backprojection alone (not a good idea in the standard imaging, where 



it leads to blurring [20||27|). We will see in the next sections that this is indeed 
the solution. Simultaneously, we will eliminate another drawback of the standard 
tomographic techniques, which is the difficulty of quantifying the confidence level of 
detection. 
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Figure 11: BP reconstruction from 3D Radon data. 700 ballistic source particles, SNR 
0.07%. Top: Cross section through source. Bottom: Isosurfaces at 55% (blue), 70% (green) 
and 80% (red) of peak value. 



3 A probabilistic discussion 

Now we investigate why it might be possible to detect geometrically small low emission 
sources by backprojection. Let us suppose that the detectors are collimated and thus 
can determine information about the incident direction of particles. Then backpro- 
jection essentially "sends" the detected particles back along their incoming straight 
lines. In particular, these backpropagated trajectories of all ballistic source particles 
will pass through the source region. The background particles or the source particles 
that have scattered before reaching the detector, will be backprojected along straight 
lines, which are different from their original zig-zag trajectories, about which we have 
no knowledge. Therefore, in the considerations that follow, we assume that the back- 
propagated trajectories of these particles are random straight lines. We also assume 
for now that the random distribution of trajectories is uniform. 

In a nutshell, our simple argument, which will be made more precise below, is that 
if the number of lines passing through a tiny region exceeds significantly the mean of 



the background, then most probably there is a source at this location (see Fig. 12). 

Let us try to quantify this and consider N random particle trajectories in a ball Bji 
of radius R. What is the probability of at least n particles out of the total N contained 



13 



Figure 12: If the local density of lines is significantly higher in the gray region, this is likely 
due to a source being present at this location. The large circle indicates the region of interest 
(ROI) rather than any physical boundary. 



in Bji to pass through a small ball Br with radius r? Let, as before, L(^^^g^ denote the 
line in M? defined by its unit normal vector a; E 5^ and a signed distance s from the 
origin, i.e. 

Thus, the lines on the plane are mapped 1-to-l to the points on the half-cylinder 
S\. X (— oc,oc), where 5^ = {a; = (cos 0, sin 6^), G [0,7r)}, and these points are 
uniformly distributed on the half-cylinder. All lines intersecting Bji correspond to 
points on S\ x (— i?, i?), and those intersecting Br are bijectively mapped to points on 
S] X (— r, r). Thus, the probability of a random line in Br passing through Br is 



P 



area(S'^ x (— r, r)) r 
area(S'| x ~ R' 



The probability that n out of the total N lines cross Br is given by the binomial 

distribution B{N,p), and is equal to ( ]p^{l —p)^~^. Since we are in the situation 

when p is fixed by the dimension of the source and N is large, the Central Limit 
Theorem applies. It follows (since the values of Np we will see will range in thousands) 
that the binomial distribution is approximated well by the Normal distribution N(fi, a^) 
with the mean fi = Np and standard deviation a = ^/Np{l — p). If it happens that 
the number of lines n crossing Br exceeds significantly (in comparison with a) the 
mean /x, the probability of this occurring just due to random reasons is very small. 
Therefore, one can be almost certain that this clustering of trajectories is the result of 
a radioactive source located at Br- Our numerical experiments of source detection by 
backprojection agree with these expectations. 



As an example, in Fig. 13 we present a typical histogram of the number of lines 
crossing a pixel in the square D — [—1,1]^ with a uniform grid of 100^ pixels. As de- 
scribed in section [2T} we generated random background and an isotropic point source. 
The particles were detected by four arrays of collimated detectors placed on the sides 
of the square. The exact incoming directions of particles were recorded. For the result 
shown below, we backprojected each particle using a line-drawing algorithm. Namely, 
each particle was back propagated along a straight line into D by adding a value of 1 
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to all pixels in D intersected by the line. After all particles were backprojected, the 
value at each pixel of D is equal to the number of lines intersecting it. On figure [13] 
we show a histogram of the pixel values after bakprojecting particles coming from a 
random background of about 10^ particles and 10^ source particles. On the far right of 
the histogram one sees the outlier contribution from the source, which deviates about 
10.3 standard deviations from the meai][^ This result complies with discussion above: 
the appearance of such a highly unlikely outlier indicates presence of a source. 




9,400 mean - 3 std mean mean +3 std 10,600 10,800 11,000 11,200 

Figure 13: Histogram of a backprojection reconstruction from data collected by four coUi- 
mated detectors, placed on the sides of the unit square. The data simulated a source at 
the point (0.3, —0.4), with 1,000,368 background particles and 1,000 source particles. The 
mean of the backprojected image is 9, 972, and the standard deviation is 99.83. The number 
of lines at the location of the source deviates 10.3 standard deviations from the mean. The 
pixels in close proximity of the source also have an elevated number of lines passing through 
them and thus are responsible for several other far to the right outliers. They would show 
as a single pick on the detection plot. 

Now we will make these statements more precise. Suppose that the area we are 
imaging is divided into Npi^. pixels. Suppose further that n lines intersect a pixel Bj, 
with radius r. We will write n = + n^, where rig and are the number of ballistic 
source and background trajectories, respectively, which cross Bj.. We will describe a 
method to determine whether there is a source located at this pixel. 

Let us choose a threshold value fc^, such that the probability of a normal variable 
to reach more than kt standard deviations above the mean is very small. E.g., kt — ^ 
or higher suffices. If an abnormally high number n of lines passing through a pixel B^ 
is detected, i.e. \i n > nt '.— ijl + ktcr, we will claim that the pixel contains a source. 
Otherwise, such claim will not be made. The probability of at least nt lines crossing 
Br due to random reasons is approximately r := 0.5 evfc^kt/V^). We also have to take 
into account the total number Npix of pixels and the possibility of such clustering of 

^There are several other outliers, deviating by 4 - 8 standard deviations. They all correspond to the 
pixels intersecting the same source. 
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lines occurring inside at least one of them. The probability of at least rit lines crossing 
at least one of the pixels due to random reasons (and thus causing a "false alarm") can 
be approximated by 

fprate=l-(l-r)^^-. (6) 

Here fp rate stands for false positive rate, which is the rate of false detections The 
above estimate assumed independence across different pixels of the events "the pixel 
is crossed by n lines". It is easy to see that this is in fact not the case: If a pixel is 
crossed by a large number of lines, there is a good chance that the surrounding pixels 
would also have a many lines intersecting. However, we claim that treating them as 
independent does not have a significant impact on the resulting estimates. To confirm 
this, we show some statistics from numerical simulations of random backgrounds at 
the end of this section. 

The true negative rate (or specificity), is given by 

tn rate = 1 — fp rate 

and represents the rate of correctly classified negative outcomes. The true negative 
rate is also the confidence probability that we have found a true source. We will write 

tn rate = confidence := (1 - r)^^'^ = [1 - 0.5 erfc(/ct/\/2)]^p^". (7) 

In particular, if the threshold value kt is set to 5 and there are 100^ pixels, the true 
negative rate is .997, which gives also this number as the confidence that the alleged 
sources are indeed true. This confidence probability depends only on the pre-set thresh- 
old value kt of the method and the number of pixels Npix- However, if in some pixel the 
number of lines deviates from the mean even further, this increases the confidence prob- 
ability that there is indeed a source there. For instance, if n > + 7cr, the confidence 
reaches, for all practical purposes, 100%. 

Let us discuss now the probability of missing a source. We assume that an apriori 
ballpark estimate of possible value of ris (the number of ballistic particles coming from 
the source) is available. We can measure Ug in the a units: 

ris — kscr. 

Here kg is a positive constant and a is the standard deviation of the distribution of the 
number of background trajectories passing through a pixeQ Now suppose that there 
is indeed a source present in a pixel B^. How can we miss it? Since our test identifies 
a source whenever n = + > + fc^a, the necessary condition for missing it is 
ni, < fi — (kg — kt)(7. The probability of this happening (i.e., the false negative rate (fn 
rate)) is given by 

fn rate = 0.5 erfc 

For instance, if kg > h -\- 3, the probability of missing this source is at most 0.0013. 
The sensitivity, or true positive rate (tp rate), is given by 

tp rate = 1 — fn rate. 

^ Recall that a is determined by the total number of detected particles, which is known, and the ratios of 
the source dimension to the dimension of the imaged region, for which we would, in practice, have a crude 
idea. 



/ kg - kt \ 
V V2 )' 
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Note that the sensitivity of the method depends on the threshold value kt as well 
as on our estimate of detected source particles. 

From the above discussion we see that in order to find sources reliably, it suffices 
to detect rig > {h + 3)cr ballistic particles from the source. 

We claim that for any given SNR 5 > and for any fc^ > it is possible to 
detect the source, provided that the total number of detected particles, 
is sufficiently large. 

Indeed, if a signal to noise ratio s is known, and if p <C 1, the number of detected 
source particles will be rig ^ sN and a ^ ^/Np{l — p). Then, we require the following 
inequality to be satisfied: 

sN > {kt + 3) ^Np{l-p). (9) 

Since the left hand side of ^ grows linearly with N and the right hand side grows only 
as \/]V, the inequality will be satisfied for a sufficiently large A^. We come up with the 
following rule of thumb: 

One can expect to detect reliably the source with high sensitivity and 
specificity, if the total particle count N at the detectors is on the order of 

iVoc f-) p{l-p), (10) 



or higher, where s is the SNR (ballistic particles from the source vs the 
total hits) and p — r / where r and R are the linear dimensions of the 
source and the whole area correspondingly. 

Let us play with some practically possible values of the parameters. A nuclear 
source with shielding could have a radius of the order of several centimeters, while a 
vehicle or a cargo container has size in the order of meters. Thus, we will assume that 

r = 0.01 ^p=10"^ (11) 



Substituting this value of p and s ^ 10~^ into 10, one expects that when N reaches 
close to 640, 000 and directional information at detectors is available, the sources should 
be reliably detectable. This squares well with the results of our previous computations 
with the backprojection methods. 

The conclusion from this probabilistic discussion is that the simple backprojection 
should reliably detect sources with low SNR, as long as the total number N of detected 
particles is high enough (see ([Tol)). E.g., for the SNR of about 0.1% the values of 
around 640,000 should be sufficient. Besides, the confidence probability of detection 
(or absence of a source) can be computed. 

The use of backprojection should also smooth out unstructured noise. In the next 
section we present numeric examples of such backprojection detections. 

It is worthwhile to notice that, as long as the threshold parameter kt (equivalently, 
the desired confidence probability) is fixed, the whole detection procedure can be done 
automatically, without necessity of a human eye assessment of the reconstruction. 

Comparison to simulation data. Recah that equation ^ and most of the fol- 
lowing discussion assumed (incorrectly) independence for different pixels of the events 
"a pixel is crossed by n lines". To con&m that the confidences we estimated under 
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this assumption are not vastly optimistic, we simulate a large number of random back- 
ground fields and compare the statistics of these samples to our theoretical estimates. 
Each sample consists of 10^ uniformly distributed random lines on a grid of 100x100 
pixels. For this setting the mean and standard deviation of the number of lines cross- 
ing a pixel are approximately fi = 10, 000 and a = 99.5. For different values of the 
threshold parameter kt ranging from 4 to 5, we record the true negative rate rt of the 
sample population, which is the ratio of samples in which each pixel had less or equal 
than rtt — jJi -\- kfcr lines intersecting. We compare this to the estimated confidences 
given by Q. The results from 50,000 samples are shown in table |3j It can be seen 
that for all values of kt the true negative rate from the sample population is slightly 
lower but very close to the estimated confidences. We conclude from this that the sim- 
plifying independence assumption made in this section is reasonable and the estimated 
probabilities agree well the true probabilities. 



h 


nt 


confidence 


n 


4 


10,398 


0.7285 


0.7141 


4.1 


10,408 


0.8134 


0.8014 


4.2 


10,418 


0.8751 


0.8658 


4.3 


10,428 


0.9182 


0.9125 


4.4 


10,438 


0.9473 


0.9428 


4.5 


10,448 


0.9666 


0.9623 


4.6 


10,458 


0.9791 


0.9760 


4.7 


10,468 


0.9871 


0.9846 


4.8 


10,478 


0.9921 


0.9907 


4.9 


10,488 


0.9952 


0.9944 


5 


10,497 


0.9971 


0.9968 



Table 2: Comparison between estimated con- 
fidences and statistics from 50,000 negative 
samples 



4 Source Detection Using Backprojection 

We illustrate the conclusions of the previous Section on a couple of typical examples 
of 2D X-ray data. The data was generated in the way described in Section |2.1[ The 
incoming trajectories of detected particles were backprojected using a line drawing 
algorithm, as was explained in Section[3| Recall that we recorded approximate locations 
of the site at which a detection occurs, since in reality detector arrays consist of a 
number of bins (detectors). 

In the first example, we simulated 639, 954 background particles and 640 particles 
from a source located at (—0.43,-0.11). All particles were backprojected and the 



source was detected with 100% confidence (see Fig 14). 

For our second example, we used only three detector arrays along three sides of a 
square. Such a "detector gate" is a more realistic configuration than detectors sur- 



rounding the whole object. Figure 15 shows a typical backprojection reconstruction. 
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Figure 14: Left: Backprojection reconstruction from X-ray data. Right: Peaks exceeding the 
local mean more than 4 standard deviations. The unit on the vertical axis is one standard 
deviation from the mean. The maximum number of standard deviations above the mean 
equals 6.3 and occurs at the location of the source. The estimated confidence is 100%. Data 
consisted of 639, 954 background particles and 640 particles coming from a source located at 
(—0.43,-0.11). The black lines represent the location of detector arrays. 



The drop of the values at one side of the square is due to detectors missing there. 
Thus, effectively, there is a non-uniform distribution of the background, which contra- 
dicts to our assumption of uniformity. However, this non-uniformity is smooth (slowly 
varying), so we can assume that locally the distribution of background trajectories is 
uniform. In order to apply our probabilistic arguments described in Section |3] we es- 
timate the local mean and standard deviation at every point. For the reconstruction 



shown on Fig, 15, this was done by analyzing the data on a 7 x 7 sub-grid surrounding 



each pixel. Then, the local means were subtracted and the resulting values plotted (see 



the bottom left picture in Fig. 15) in the units that are equal to the local standard 



deviation a. Finally, the values below the threshold fc^ = 4 were erased (bottom right 



in Fig. 15). The result clearly shows the detected source. The confidence probability 



was calculated as 

confidence = (1 — 0.5 * erfc^ksource/V^))^^^^ , 

where kgource is the height of the highest peak. The calculated confidence level for this 
example was 99.99%. Notice that if we used the threshold value fc^ = 4 instead of the 
actual peak's height, we would have significantly underestimated the confidence. 



5 Compton Cameras 



Compton cameras, also called Compton scatter cameras or electronically collimated 
cameras, have received a lot of attention since they were first proposed in [41j and a 
prototype was developed in ISGIIST). Compton cameras were first suggested for use in 



nuclear medicine (SPECT), but they also have applications in astrophysics 14,38,44 
monitoring nuclear power plants |33,34, and other areas. Compton cameras offer 
several advantages over the conventional mechanically collimated (Anger) 7-cameras. 
The most significant of these is the dramatic increase of sensitivity - Compton cameras 
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Figure 15: Backprojection reconstruction from X-ray data (top). The unit on the vertical 
axis in the bottom pictures is one (local) standard deviation from the (local) mean. Bottom 
left picture shows the result of backprojection after subtraction of the local mean. Peaks 
that exceed the local mean by more than 4 (local) standard deviations are shown at bottom 
right. Source located at (0.203,0.101) is detected with confidence 99.99%. Data consisted 
of 639,417 background and 646 source particles. The black lines represent the location of 
detector arrays. 



are reported to count orders of magnitude more photons than Anger cameras 



Another advantage is the flexibility of geometrical design [33|[34], which even allows 
for hand- held devices [T7l[23] . 

A Compton camera consists of two detectors, which, in the standard setup, are 
planar and are placed one behind the other (see Fig. 16). A photon incident on the 
camera undergoes Compton scattering in the first detector, which records the location 
xi and the energy of interaction Ei. After the scattering, the photon is absorbed in 
the second detector, where the position X2 and energy of absorption E2 are measured. 
From the knowledge of Ei and E2, the scattering angle '0 can be determined by the 
formula (e.g. [7||4T]) 

mc^Ei , 

Here m is the mass of an electron and c is the speed of light. The direction into which 
a photon scatters is given by —/3 := \. From the knowledge of /3 and the 

|X2 - Xi| 

scattering angle i(j we conclude that the photon originated from the surface of the cone 
with central axis /3, vertex xi and opening angle 2t/j. Therefore, although the exact 



20 




Figure 16: Schematic representation of a Compton camera. 



incoming direction of the detected particle is not available, one knows a cone of such 
possible directions. The goal of Compton camera imaging is to recover the distribution 
of the radiation sources from this data. 

The description above addresses the Compton 7-cameras only. However, neutron 
detectors are currently being developed (25j|40] that (employing a different physics 
principle) provide the similar cone information. We thus will not make distinction 
between these and call them all "Compton type cameras." 

We will now address briefly the mathematics behind the Compton type measure- 
ments. 



5.1 The Cone Transform 

Suppose that we want to image the distribution /(y) of radioactivity sources inside 
a certain object in R^. From now on, we will assume that /(y) is supported on one 
side of the Compton camera and that supp/ does not intersect the camera. Under the 
assumption of a sufficiently long observation time (and thus a large number of particles 
detected), the Compton camera provides us with the projections [l]|7||43!|, which we 
denote by C/(x,/3,'0), i.e. the integrals of /(y) over cones parameterized by a vertex 
X lying on the detector, central axis vector j3 from the unit sphere 5^, and half-angle 



7/^ G [0,7r] (see Fig, 17): 



Cf{yi,f3,i;) = K{i;) / / /(x + ra((/)))r sin^ dr a-^ = cos^, a G 5^ (13) 
Jo Jo 

Here K{t/;) is the Klein-Nishina distribution of scattering angle^ The function C/(x, /3, tp) 
provides the expectation of the number of hits, per unit time, from the particles moving 
along the cone. In particular, when the particle count is low, these values are not well 
determined by the data. This issue was addressed in the preceding Sections. An im- 
portant observation is that the problem of inverting the cone transform C/(x,/3,'0) is 
even more overdetermined than the 3D X-ray transform. Indeed, /(y) is a function of 
three variables, while Cf('K^/3^il;) depends on five parameters (here x is restricted to a 



^The Klein-Nishina factor is a known function and thus can be easily accounted for. We therefore do not 
include it in the further analysis. 
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Figure 17: A Compton cone with apex x, central axis f3 and half-angle ip. 



2D detector surface). Most authors consider restricted versions of the cone transform, 



reducing the number of parameters from five to three |lj|7,29 . The paper |39| contains 
a closed form inversion formula, which uses four of the parameter^ We, however, have 
to use the full 5-dimensional dataset, since, as it was mentioned before, any attempt 
to restrict the data would wipe out the signal. 

5.2 Reconstruction Techniques in Compton Camera Imag- 
ing 

We provide here a brief survey of some known reconstruction methods in Compton 
cameras imaging. Many of them reduce cone projections to Radon projections, i.e. 
integrals of / over planes. 

Algebraic Reconstruction Techniques, as elsewhere in tomography, are widely used 
in Compton cameras imaging. This is partially due to the fact that analytic inver- 
sion formulas were not available until the second half of 1990's. Maximum likelihood 
methods were developed in [3||8|[l^. A maximum likelihood method was also used 
in imaging with the COMPTEL telescope |38|. In 1 35 a matrix inversion technique 
was utilized to solve the problem for a specific detector geometry. Fast backprojection 
algorithms were developed in f32^,'45|. 

A spherical harmonics expansion solution was proposed in |1 . For every point 
x on the detector array, and a fixed in advance half-angle '0 the authors relate the 
conical projection C/(x, /3,'0) to the Radon projection (integral) of / over the plane 
passing through x and perpendicular to /3. A fast algorithm for computing the spherical 
harmonic series was also developed in this paper. Several other authors provided 



spherical harmonics solutions, however their model for Compton data differs from (13), 
e.g. (I9l|3l}|42]. 

Probably the first closed form analytic inversion formula was obtained in |7|. The 
authors considered only cones perpendicular to the detector plane, i.e. cones with 
central axis /3 = z. Thus, the Compton data depends only on the vertex x := (x, y) and 
the opening angle tp. Let us denote g{x, y, ip) \— C f{'K, z, ip) and let F{u, v, •), G{u, •) 
be the two-dimensional (x, y)-Fourier transforms of /(x, z) and g(x, ip) respectively. 



^Recently [26 , the complete set of data was used for reconstructions. 
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Then the fohowing relation holds: 

--G{u, V, t) 



(14) 



where ^ = zVu'^ + v'^ and t = tan^. Here T-Lq is the zero-order Hankel transform 

roo 

T-Log{p) = 27r / g{r)rJo{27vrp) dr, 
Jo 
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and Jo (7) = (27r)-i / e^^^"""^ # is the zero-order Bessel function. Later, in 
^0 

more properties of the same restricted cone transform were given. 

In [26] the authors use yet another mathematical formulation for the forward prob- 
lem, which accounts for the efficiency of the detector at different incident angles of 
particles. They extended the approach of |7| to show that the radioactivity distri- 
bution can be reconstructed from any family of cones that have central axes /3 such 
that the angle between (3 and the z-axis is fixed. Then averaging of the solutions for 
different values of this angle was employed, thus making use of all available data and 
reducing noise. 

The following useful relation between the cone and Radon transforms was found 
in (391: 



1 roo 1 PTT 

HRf{f3,^'f3) := - / Rf{f3,t)h{^-f3-t) dt = / C/(x, ^, ^)/i(cos^) # (15) 

^ J-oo ^ Jo 

Here H denotes the Hilbert transform and R is the Radon transform in three dimen- 
sions. Note that four out of five variables are used and that the i7i?/(/3,x • /3) is 
constant on the line Lc := {x E detector plane|x-/3 = c}, for any fixed constant c. The 
last observation shows that it is possible to reconstruct / if the vertices of cones x are 
restricted to a curve. This also allows for averaging of the solution on the lines Lc- In 



Section 6.1 we will consider a two-dimensional analog of this formula. 



6 Compton Cameras in Two Dimensions 

In the two dimensional setting we assume that the detectors x lie on a line, which we 
call the detector line. A cone in two dimensions is defined by a point x that serves as 
its vertex, a central axis unit vector /3, and a half-angle '0. Such cones simply consist 
of two rays with a common vertex. More precisely, let /3 = (cos /3, sin /3) G S^-*^, /3 G [0, tt] 
be the central axis of a cone with opening half-angle G [0, tt]. The two half-lines that 
constitute the cone pass through the detector point x and are given by 

{yGR2:y = ra„r>0},z = l,2, 

where cxi = (cos(/3 — '0), sin(/3 — ip)) and CX2 = (cos(/3 + '0), sin(/3 + '0)), see Fig. 18 
As we have agreed before, we drop the known factor K(t(;). Thus, in what follows we 
will assume uniform distribution of scattering angles, that is K{t/j) = 1. Then the two 
dimensional cone transform Cf{'K,/3,il;) is the projection of /(y) onto these lines: 

roo 

C/(x, (3,ij)^ / (/(x + rai) + /(x + ras)) dr = 2?/(x, ai) + I?/(x, 02). (16) 
Jo 
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detector line 3 

X 



Figure 18: Two-dimensional cone with apex x, central axis /3 and half-angle ip. 



Here I^/(x, a) = /(x + ra) dr denotes the divergent beam (or fanbeam) transform 
of /. 

For convenience, we will sometimes express Cf and Vf as functions of the angle /3 
instead of the unit vector /3. Thus, we can write (flGl) as 



C/(x, /3, i,) = Vfix, + Vf{x, (3 + ^P). (17) 

As in 3D, the problem of inverting the two-dimensional transform is overdeter- 
mined: Cf{'K,/3,i(;) depends on three parameters, while /(y) only on two. 

The authors are aware of few papers devoted to Compton cameras in two dimen- 
sions. In |2,13| subsets of the cone projections of f{x) are represented as line integrals 
of the function / added with its mirrored shear transformation. Namely, let us intro- 
duce the variables ki = cot(/3 + '0), k2 = cot(/3 — '0). The cone transform can be written 
as 

roo roo 

Cf{x,(3,j/j)= / f{x + kiz,z) dz^ / f{x^k2Z,z) dz, (18) 
Jo Jo 

where the x— axis coincides with the detector line and the z— axis is perpendicular to 
the detectors. Let us fix the sum ki-^ k2 — K and define the function 

r ( ^ if{x,z), if z > 

\j[x — Kz, — z), ir z < 0. 

Then, Cf{x,l3,ijj) = J^^fxix + kiz^z). Thus, for each fixed value of the cone 
projections are represented as line integrals of a certain function, so all that remains 
is to invert the Radon transform. The reconstructed image will contain the function / 
above the detector line and a mirrored shear transformation of / below the detectors. 
Obviously, by fixing the parameter K, not all available data is used. The case when K = 
corresponds to using only those cones with central axes orthogonal to the detector 
line. Since we need to use all over deter mined data, this approach is inconvenient. 
Besides, as our previous consideration of tomographic methods shows, in our specific 
situation, backprojection methods, rather than attempting the "full" reconstruction 
should be used. 
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6.1 Three Inversion Methods 



6.1.1 Method A 

The following relation, proven in makes use of all parameters and is an analog of 
the three dimensional result presented in [39) . 

Theorem 4. Let the closed unit ball lie on one side of and away from the detector 
line, and let f{x) E Hq^'^{B'^). Denote h{t) = l/t. Then 

-1 roo -1 riT 

Hnf{^-f3,f3) = - / nf{f3,t)h{^-f3-t) dt = — / C/(x,/3,^)/i(cos^) d^/j (20) 
^ J-oo ^ Jo 

Here H is the Hilbert transform in the first (scalar) variable and the integrals are 
understood in the principal value sense. 

Here is the unit disk in the plane and Hq is the standard notation for Sobolev 
spaces (e.g., |9|). 

Corollary 5. Theorem^ provides an inversion formula for the cone transform. 



Indeed, computing the integral in the right hand side of (20), one recovers HTZf ^ 
where TZ is the 2D Radon transform, and H is the Hilbert transform with respect 
to the linear variable. Then, the filtered backprojection formula for inversion of the 
two-dimensional Radon transform, 

/(-) = i^* {^i^f) (21) 



implies that after differentiating the right hand side of (20) with respect to the linear 
variable and backprojecting, one recovers the function /. 



Remark 6. In the left hand side of (20), one sees TZf{s, j3), where s — i^- j3. Thus, if 
one infinite detector array is used and (3 is perpendicular to the array, we would only 
know 7lf{s,/3) for a single value of s. If the detector array is finite, as is the case in 
all implementations, we would miss a much bigger chunk of data. One possible way to 
obtain the complete Radon data in {s,/3) G [—1,1] x is, for example, to use three 
finite size detector arrays placed along the sides of a square containing the object. 

6.1.2 Method B 

Another, simpler way of obtaining Radon data from cone data is presented below. Let 
us denote by u(k, a) the integral of / along the ray starting at the point x in direction 
of the unit vector a (i.e. the fanbeam projection of /): 



POO 

u(K,a.)=Vf(K,a.):= / f(K + rcx)dr. 

Jo 



(22) 



Let us fix a detector position x and a direction (Xq and try to recover the fanbeam data 
ix(x, ao). For any G [— tt, tt] we can write C/(x, /3, = ix(x, /3 — +i/(x, /3 + 
Thus, the following relation holds 

i/(x, ao) = C/(x, ao + '0, |'0|) - ix(x, ao + 2^) 
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Figure 19: In order to determine the integral of / along the line Iq^ add the integral over 
the cone consisting of Iq and li to the integral over the cone determined by Iq and l2- Then 
subtract the integral over the cone determined by li and I2. 



It is easily seen now (Fig. 19) that for any two half-angles '0i,'02 ^ [—^5^] the 
fanbeam data i^(x, 0:0) can be obtained from cone data as follows: 

i/(x,Q:o) 

^[C/(x,ao + ^i,|^i|) + 

C/(x, ao + ^2, 1^2!) - C/(x, ao + ^1 + ^2, 1^1 - ^2!)] 

(23) 

The angles ^1 and i(j2 were arbitrarily chosen, so one could average on i/ji and '02 in 
order to use all available data. This averaging also helps reducing the effects of the 
background noise. 

If averaging over half-angles '0i and '02 is used, this method is computationally 
expensive. We propose a third reconstruction method, which is fast and makes use of 
all available data. 



6.1.3 Method C 

Let us recall that the source intensity function /(y) is compactly supported and supp/ 
lies on one side and away from the detector array. Let us fix a detector location a. Then 
the function i^(a, a) (see ([22])) is supported in (0, tt), so it can be extended periodically 
in a from [0, 27r] to (— oc, 00). Let us define 

/7T—a 
C/(a,a + V, IV'I) # (24) 
-a 

We can write, 

/7T — Q. PTT — Q. 

i/(a, a) (i'0 + / i/(a, a + 2'0) di(j} = 

.2. (25) 
u{a, a) + (27r) ^ / u{a, a) da 
Jo 

In the above equation, the periodicity of u{a, a) with respect to a was used. Then, 

u{a, a) = u{a, a) + ^ jj^,^^ /(a + rw) drdu; 

= P/(a,a) + (47r)-i(7e#7e/)(a). ^ ^ 
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Here, as before, VJ^ denotes the backprojection operator 

T^*9{y) = / ^(y • ^,^) dio. 

j\u\=i 

The second term in the right-hand side of ([26]) can be written as 

P[(47r)-i7e#7e/(a)5(a-y)], 
where 5 is the Dirac's delta-function. Therefore, we come to the fohowing conclusion: 

Theorem 7. When the inverse fan-beam transform is applied to u, the result equals 
the sum of the function /(y) and a distribution supported on the detector array only. 
Hence, functions f supported away from the detector arrays are recovered correctly. 

In essence, this method approximates an integral of a function on a given line by 
averaging all its cone integrals which have one side lying on the line. 

7 2D Compton data backprojection examples 

In order to generate Compton camera data, we first generated X-ray data in the manner 
described in Section [2T| Recall that the X-ray data for each particle consists of a bin 
on a detector array and exact incoming direction a. Then, for every particle we chose a 
scattering angle ^ drawn from a uniform distribution in (—a, n — a). The central axis 
of the scattering cone was computed as the sum j3 — a^ifj. The information provided 
by the Compton camera is three dimensional and consists of the bin on the detector 
array, the central axis of the cone /3 and the absolute value of the scattering angle '0. 

For the reconstructions from Compton camera data presented in this section, we 
used Method C, described in Section [6.1. 3[ to convert the cone data into X-ray data. 



Method C was chosen as the least computationally expensive and also using all cone 
data. Note that for the particular type of data we are handling, the central axis and 
scattering angle of cones are continuous random variables. Thus, for a finite number 
detected particles, one cannot average on cones with a common side, because such 
cones are a rare occurrence. Therefore, the method is equivalent to simply treating a 
particle coming from a cone (x, /3, ^) as two separate particles detected at x and having 
incoming directions 13 — and (3 -\- tp. After converting cone data into X-ray data in 
this way, we employed the usual reconstruction procedures for X-ray data. 

It is also important to make clear that in our experiments we assume that exact cone 
axes and scattering angles are known. The effects of uncertainty of these quantities 
are to be investigated in the future. 

First, we used filtered backprojection to reconstruct the source distribution from 
data provided by four Compton detector arrays placed on the sides of the square [—1,1]. 
The signal to noise ratio was fixed at 0.1% and we varied the bulk number of detected 
particles. As for the case of X-ray data, we defined a detection successful provided that 
highest peak of the reconstructed image occurred at the location of the source. The 
number of successful FBP detections increased slowly as the bulk number of detected 
particles went from 600,000 to 1,000,000. The results are reported in Table [7| An 
example of a successful detection using Method C combined with FBP is shown in 
Fig. I20| 
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method \ bkgd particles 


600, 000 


700, 000 


800, 000 


900, 000 


1,000,000 


Method C + FBP 
BP + local grids 


8/20 
17/20 


14/20 
17/20 


12/20 
20/20 


18/20 
20/20 


17/20 
19/20 



Table 3: Number of successful detections from cone data of a source located at 
(0.311,-0.433), with a fixed SNR 0.1% and varying number of detected particles, using 
different reconstruction methods. For each level of background particles 20 sets of random 
data were generated. In each case we applied Method C combined with FBP and backprojec- 
tion followed by local grids of size 9x9. In all but one experiment, if an FBP detection was 
successful, so was the corresponding backprojection detection. The one example in which 
FBP detection occurred but backprojection failed was at the level of 700, 000 background 
particles. 




Figure 20: Filtered backprojection reconstruction from Compton data (left). The number 
of background particles is 900020 versus 900 source particles. The black lines represent 
detector arrays. The right picture shows only the peaks which deviate more than 4.3 standard 
deviations from the mean. The correct source location is found. Pixels close to the detectors 
are not shown, as the values of the reconstructed image there are quite high, see Theorem [7| 



Next, we use Method C for converting cone data into X-ray data, followed by 
backprojection based on the line drawing algorithm described above. Thus, for every 
detected particle, two particles are backprojected along the sides of the cone. Therefore, 
the value of the resulting image at every pixel equals the number of cones intersecting 
it. OnFig.[2l]^top) we show a typical backprojection reconstruction from four Compton 
detector arrays from exactly the same data used for the previous example (Fig. 20). 
The background is not uniformly distributed, so we employ the local grids approach 
introduced in Section |4| The size of the local grids was chosen to be 9 x 9 pixels. On 
Fig. [2T|bottom) a plot of the number of local standard deviations away from the local 
mean is shown. The value at the source location deviates 4.93 standard deviations 
from the mean, and the estimated confidence is 99.59%. 

In Table [7| we compared FBP reconstructions with backprojection reconstructions. 
For all backprojection results we used local grids of size 9x9 pixels on which the local 
means and standard deviations were estimated. In order to compare backprojection to 
FBP, we deemed a detection successful, if the number of local standard deviations was 
highest at the location of the source. Our results show that backprojection outperforms 
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Figure 21: Unprocessed backprojection reconstruction from Compton data (top). The num- 
ber of background particles is 900020 versus 900 source particles. The bottom left shows 
deviations from the (local) mean, measured in the units of (local) standard deviations. The 
size of local grids is 9 x 9 pixels. Only the peaks that exceed the local mean more than 4.3 
(local) standard deviations are shown on bottom right. Source located at (0.311, —0.433) is 
detected with confidence 99.59%. 



FBP, as in the case of X-ray data. The results are summarized in Table [7| 

As we have already seen, detection based on backprojection allows for estimation of 
confidence. Similarly to detection by X-ray backprojection, one could set a threshold 
for detection at 4.1 or 4.3 local standard deviations above the local mean. According 
to our estimate for the confidence ([7|), such thresholds would result in confidences at 
least 81% or 91.8% respectively. In Table [t] the number of successful detections for 
these two values of the threshold parameter are given. We used the same simulation 
data as for Table [71 



threshold \ bkgd particles 


600, 000 


700, 000 


800, 000 


900, 000 


1,000,000 


4.1 stds above mean 
4.3 stds above mean 


10/20 
7/20 


9/20 
9/20 


14/20 
11/20 


18/20 
14/20 


18/20 
16/20 



Table 4: Number of successful detections from cone data with a fixed SNR 0.1% and varying 
number of detected particles. Backprojection was used and local means and standard devia- 
tions were estimated on local grids of size 9x9. Two diflFerent threshold values for detection 
were set, resulting in confidences at least 81% and 91.8% respectively. 

In comparison with X-ray data, a larger number of background particles is needed 
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for detection with Compton cameras. This reflects the somewhat lower quality of 
Compton data in comparison with the data from collimated detectors. 

A couple of examples of reconstructions from three Compton cameras, which form 
a "gate" of detectors, are provided below. 



The source in Fig. [22] is close to the side of the square where there is no detector, 
which makes it harder to observe. The number of ballistic source particles was 1050 
versus 999, 925 background particles. The local mean and standard deviation were 
estimated on 9 x 9 local grids. The maximum number of local standard deviations was 
5.56 and occurred at the location of the source. The estimated confidence of detection 
is 99. 





i 




Figure 22: Backprojection reconstruction from Compton data (top). Detected source par- 
ticles - 1050, background particles - 999,925, source location - (—0.503,0.11). Number of 
(local) standard deviations above the (local) mean at each pixel is shown on bottom left. 
The size of local grids is Peaks that deviate more than 4.3 (local) standard deviations from 
the mean are shown on bottom right. The source is detected with confidence 99.98%. 



Finally, we present an example of a successful detection of several sources (Fig. 23). 
Three sources were placed inside the imaged area, and three Compton-type detector 
arrays were used to detect particles. The total number of detected ballistic source par- 
ticles was 3, 229, each of the sources having roughly the same contribution. Additional 
999, 325 random background particles were detected. We have set up the threshold to 
kt = 4.3a, which gives confidence of detection at least 91.97%. Note that this confi- 
dence probability is estimated based on the threshold, and not on the actual number 
of deviations at each peak, as in the previous examples. Thus, individual confidence 



30 



probabilities for the three sources discovered are even higher. The local grids we used 
for estimating the local means and standard deviations had dimensions 9x9. 




Figure 23: Backprojection reconstruction from Compton data (top). Sources located at 
(0.203,-0.101), (-0.301,0.43) and (0.41,-0.61) emitted 1108, 981 and 1140 particles, cor- 
respondingly. The number of detected background particles was 999, 325. Number of (local) 
standard deviations above the (local) mean at each pixel is shown on bottom left. Peaks 
that deviate more than 4.3 (local) standard deviations from the mean is shown on bottom 
right. Confidence of detection based on the threshold kt = 4.3 is 91.97%. 



8 Final remarks and conclusions 

Here are the main conclusions that we have reached: 

1. Direction sensitive detectors are needed to be able to uncover existence of a low 
emission source in the presence of dominating background noise. Collimated 
cameras, although providing the most valuable direction information, are not 
applicable, since collimation would decimate the already low signal. Compton 
type cameras can be used instead. 

2. Although standard integral- geometric models of emission tomography (X-ray and 
Radon transforms and their attenuated versions) are not applicable to the situa- 
tion of a very weak source, they work to some extent in the detection problem. 
However, they fail way before reaching the low values of SNR needed in applica- 
tions. 
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3. In fact, only the backprojection part of the tomographic reconstruction does the 
detection work, and the filtering part worsens the results. 

4. The (unknown) attenuation might lower the SNR, but otherwise is irrelevant for 
the validity of the backprojection method. 

5. The conclusions above hold both for the collimated and Compton type cameras. 

6. The backprojection procedure can be justified by a simple probabilistic consider- 
ation, which also provides confidence probabilities of detection. 

7. A simple algorithm is provided that recovers the source from Compton data at 
SNRs on the order 0.1% and beyond. It also provides confidence probabilities. In 
theoretical considerations and numerical tests, the algorithm demonstrates high 
sensitivity and specificity. The procedure can be made completely automatic, 
without the need of a human eye assessment of the image. 

8. The algorithm is based upon the assumption of an uniform random background. 
However, the non-uniformity that arises when the object is only partially sur- 
rounded by detector arrays, can be handled successfully by using local grids. 

9. As it is done in SPECT, one can try to use Bayesian methods to approach the 
same detection problem. This was done in 2D case in |46 . 

Certainly, there remain quite a few questions that need to be resolved. Some of them 
will be discussed in the next publication. In particular, 3i? Compton backprojection 
algorithms will be developed and studied analytically and numerically; more complex 
cargo structures will be modeled and tested, which will involve non-uniform attenuation 
and scattering; local grid version of the algorithm will be tested on more complex non- 
uniform random backgrounds than the ones arising due to the partial view; dependence 
of the results on the precision of Compton type cameras will also be addressed. We 
assumed that the exact incoming direction of particles is detected. In practice this is 
usually not the case, as only a probability distribution of the directions may be known. 
The effect of directional uncertainty is important from practical standpoint and we 
plan to consider it in the future. 
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